Notas de clase: Sistemas Lineales

Table of Contents
Notas de clase: Sistemas Lineales Clase 1: Presentación del curso Clase 2: Repaso de ecuaciones diferenciales ordinarias Definiciones Ejemplos 1) La ecuación diferencial como definición de un campo vectorial 2) Solución de ecuación homogénea 3) Solución de ecuación no homogénea Clase 3: Ecuaciones en diferencias Definiciones Ecuación en diferencias a partir de una ecuación diferencial: Ejemplo: Ecuación en diferencias hacia atrás de orden 2 Comparación entre solución discretizada y la solución de la ecuación diferencial discretizada T implícito vs T explícito Solución analítica de una ecuación en diferencias Primer ejemplo de solución analítica Segundo ejemplo de solución analítica: lambda imaginario Tercer ejemplo de solución analítica: Ceros que no tienen sentido Clase 4: Transformada de Laplace Definición: Transformada de Laplace Transformadas básicas transformada de Laplace Propiedades básicas transformada de Laplace Clases 5 y 6: Transformada Z Concepto de función muestreada Transformada Z Transformadas básicas Propiedades operacionales Transformada Z inversa Solución de ecuaciones en diferencias aplicando transformada Z en MatLab Solución mediante fracciones parciales Solución mediante la serie infinita de potencias Clases 7 y 8: Introducción a las funciones de transferencia Definición de la función de transferencia Conceptos claves Clases 9 y 10: Funciones de transferencia Retardos y orden relativo Sistemas de fase no mínima El efecto de los ceros de la función de transferencia Comportamiento de un sistema de fase no mínima Reducción de orden Cancelación de polos y ceros de fase mínima Eliminación de polos insignificantes Ejemplos de reducción de orden Caso especial de reducción de orden Función de ponderación (Secuencia de ponderación) Aproximación de la respuesta temporal mediante la secuencia de ponderación Discretización de la función de transferencia Transformada z modificada Ejemplos de discretización Clase 11: Introducción al espacio de estado Definiciones Ecuaciones de estado Ecuación en el espacio de estado continuo Ecuación en el espacio de estado discreto Obtención de la ecuación de estado a partir de la ecuación diferencial Para sistemas continuos Para sistemas discretos Solución de la ecuación lineal de estado continua Solución mediante valores y vectores propios Solución mediante serie de potencias (matriz exponencial o de transición del estado) Solución mediante transformada de Laplace Ejemplo Solución de la ecuación lineal de estado discreta Solución iterativa Solución mediante transformada Z Ejemplo Transformaciones lineales y formas canónicas Ejemplo Función de transferencia MIMO Ejemplos Discretización en el espacio de estado Modelado de sistemas no lineales Características de los sistemas no lineales Principio de superposición Puntos de equlibrio Bifurcación Caos Respuesta a entradas sinusoidales Comportamiento asintótico Funciones adicionales para el curso

Clase 1: Presentación del curso

Clase 2: Repaso de ecuaciones diferenciales ordinarias

clear,clc

Definiciones

  1. Ecuación diferencial: Ecuación que contiene derivadas
  2. Ecuación diferenciales ordinarias: Aplicación de las ecuaciones diferenciales a sistemas donde los parámetros son concentrados (todas las propiedades del sistema se resumen a un punto, es decir, se minimizan las variables independientes necesarias). Ejemplo:
  3. Ecuaciones diferenciales parciales: Aplicación a sistemas con parámetros distribuidos (las propiedades del sistema requieren más de una variable para expresarse)
  4. Ecuaciones lineales: Se dice que una ecuación diferencial de n-ésimo orden es lineal si es lineal en . Dicho de otra forma, la ecuación diferencial de n-ésimo orden es lineal si tiene la forma .
  5. Ecuaciones lineales con coeficientes constantes: Se les conoce como LTI (Linear-Time-Invariant). Los coeficientes que acompañan a la variable dependiente no cambian con el tiempo, por tanto, son de la forma Existen métodos analíticos para resolver estas ecuaciones. El término tiene significado como la entrada del sistema.
  6. Problema de valor inicial: Es un problema de la forma
  7. Solución general de un problema de valor inicial: Se calcula la solución general como la suma de la solución de la ecuación complementaria homogénea () y una solución particular de la ecuación no homogénea (). De forma tal, la solución general viene dada por .
  8. Solución de la ecuación homogénea (): Se obtiene el polinomio característico y se resuelve para obtener los valores característicos. Debe lidiarse con los casos de valores característicos reales diferentes, complejos diferentes y raíces múltiples. El polinomio característico de es de la forma .
  9. Solución de la ecuación no homogénea: Se puede resolver aplicando el método de coeficientes indeterminados (método especial pero simple) o el método de variación de parámetros (método general pero complejo). Para ambos métodos consultar el libro de Dennis G. Zill & Warren S. Wright.

Ejemplos

1) La ecuación diferencial como definición de un campo vectorial

t=-3:0.5:12; y=-3:0.5:12; % define grid of values in t and y directions
[T,Y]=meshgrid(t,y); % creates 2d matrices
dT=ones(size(T)); % dt=1 for all points
dY=27-3*Y; % dy=(t^2-y)*dt: this is the ODE
N=sqrt(dT.^2+dY.^2); % magnitude of arrows
dT=dT./N; % normalize arrows to get all same length
dY=dY./N;
figure(1)
clf
quiver(T,Y,dT,dY) % draw arrows (t,y) --> (t+dt,y+dy)
axis equal
axis([-4 13 -4 13]) % repeat adjustments...
grid on
hold on
C1=8;
y=9-C1*exp(-3*t);
plot(t,y,'r','LineWidth',2);
hold off
xlabel('t')
ylabel('y(t)')

2) Solución de ecuación homogénea

s1.png
Comprobación de la solución en MatLab
syms y(x)
ode=diff(y,6)+4*diff(y,5)+7*diff(y,4)+6*diff(y,3)+2*diff(y,2)==0
ode(x) = 
ySol(x) = dsolve(ode)
ySol(x) = 

3) Solución de ecuación no homogénea

s2.png
Comprobación de la solución en MatLab
ode2=diff(y,2)+4*diff(y)+4*y==x+2*exp(-2*x)
ode2(x) = 
ySol2(x) = simplify(dsolve(ode2))
ySol2(x) = 
Solución con condiciones iniciales
y_0=2; y_1=3;
conds = [y(0)==y_0; y(1)==y_1];
ySol3(x)=dsolve(ode2,conds)
ySol3(x) = 

Clase 3: Ecuaciones en diferencias

clear;clc

Definiciones

  1. Ecuación en diferencias: Es una ecuación que contiene diferencias (Δ).
  2. Diferencias hacia adelante: ,
  3. Diferencias hacia atrás: ,
  4. Forma canónica de una ecuación en diferencias hacia adelante: Se obtiene luego de simplificar Δ y reemplazar T,
  5. Forma canónica de una ecuación en diferencias hacia atrás: Se obtiene luego de simplificar Δ y reemplazar T,

Ecuación en diferencias a partir de una ecuación diferencial:

Supongamos que deseamos discretizar la ecuación , para ello debemos notar que (diferencias hacia atrás) y por ende podemos reescribir la ecuación como , de forma tal, el problema original se convierte en

Ejemplo: Ecuación en diferencias hacia atrás de orden 2

Resolver la ecuación en diferencias dada por
Solucionamos aprovechando la relación de recurrencia:
Nota: Es posible utilizar más de un método de discretización para aproximar las derivada, por ejemplo el método de Runge-Kutta4 o el método de Heun.

Comparación entre solución discretizada y la solución de la ecuación diferencial discretizada

Supongamos que tenemos la ecuación , al discretizar con T llegamos a la expresión . Para resolver esta ecuación en el intervalo podemos implementar el siguiente código
T=0.21;
a=1;
y_0=-4;
time=0:T:10;
 
%%%% Solución continua
syms y(t)
ode = diff(y)+a*y-exp(-t)==0
ode(t) = 
cond = y(0)==y_0;
ySol(t)=dsolve(ode,cond)
ySol(t) = 
 
%%%% Solución discreta
yDis=t;
yCont=yDis;
yDis(1)=y_0;
for k=[1:length(time)-1]+1
yDis(k)=(1-a*T)*yDis(k-1)+T*exp(-T*(k-1));
end
 
%%%%%%% Comparación entre soluciones
figure(1)
clf
fplot(ySol,[0 10])
grid on
hold on
stem(time,vpa(yDis))
xlabel('t')
ylabel('y')
legend({'Continuous','Discrete'})

T implícito vs T explícito

Una vez se ha obtenido la ecuación en diferencias a partir de la ecuación diferencial es posible omitir el instante en el tiempo y concentrarse solo en el número de muestreo actual (T implícito). Para graficar adecuadamente es necesario recordar que .
s3.png

Solución analítica de una ecuación en diferencias

De la misma forma que se pueden solucionar las ecuaciones diferenciales LTI, es posible plantear un problema de autovalores para las ecuaciones en diferencias. Sea tal solución , reemplazamos en la ecuación en diferencias para obtener , ahora, introduciendo el cambio de variable tenemos que , lo cual se cumplirá para . Resta notar entonces que para encontrar las soluciones de la ecuación.

Primer ejemplo de solución analítica

s4.png
s5.png

Segundo ejemplo de solución analítica: lambda imaginario

s6.png

Tercer ejemplo de solución analítica: Ceros que no tienen sentido

s7.png

Clase 4: Transformada de Laplace

Definición: Transformada de Laplace

s8.png

Transformadas básicas transformada de Laplace

s9.png

Propiedades básicas transformada de Laplace

s10.png

Clases 5 y 6: Transformada Z

Concepto de función muestreada

A partir de una función cualquiera , definiremos la función muestreada con intervalo de muestreo T como se sigue
Debe ser claro que

Transformada Z

Aplicando la transformada de Laplace a una función muestreada obtenemos que
, ahora, haciendo el cambio de variable , llegamos a que
, adicionalmente, sabiendo que T actúa como una variable muda se sigue que podemos asumir y considerar

Transformadas básicas

Propiedades operacionales

Transformada Z inversa

Notando que la transformada Z de una función es una fracción propia, es posible calcular la transformada Z inversa mediante procedimientos:
  1. Método de la serie infinita de potencias: Se procede mediante división larga. Los coeficientes de la división corresponden con los valores de la función en cada instante de muestreo.
  2. Desarrollo en fracciones parciales: Se debe garantizar que cada fracción parcial sea una función propia. De forma tal, buscamos reconocer estructuras a partir de las cuales aplicar la tabla de transformadas básicas y propiedades operacionales.
  3. Integral de inversión: Es posible calcular la transformada inversa por definición, sin embargo, será necesario resolver una integral compleja por el teorema del residuo de Cauchy. Este método no se utiliza en clase.

Solución de ecuaciones en diferencias aplicando transformada Z en MatLab

Supongamos que queremos resolver la siguiente ecuación en diferencias sujeta a condiciones iniciales

Solución mediante fracciones parciales

Solución analítica
clear;clc
%y(k+3)+0.8*y(k+2)==0 % y(0)=1, y(1)=0, and y(2)=1.
syms y(k) z %se declara la función simbólica y(k) y la variable simbólica z
assume(k>=0 & in(k,'integer')) % se establece que k es entero mayor o igual a 0
f = y(k+3)+0.8*y(k+2)==0 % se declara la ecuación en diferencias
f = 
fZT = ztrans(f,k,z) %transformada z de la ecuación en diferencias, la variable k pasa a z
fZT = 
syms pZT %variable para resolver la ecuación
fZT = subs(fZT,ztrans(y(k),k,z),pZT) % se sustituyen las expresiones de Y(z) por la variable simbólica pZT
fZT = 
pZT = solve(fZT,pZT) % se soluciona la ecuación para pZT=Y(Z)
pZT = 
pZT = subs(pZT,[y(0) y(1) y(2)],[1 0 1]) % se sustituyen los valores de las condiciones iniciales para obtener la fracción Y(z)
pZT = 
pZT=partfrac(pZT) % esta línea de código ayuda a evidenciar la descomposición en fracciones parciales
pZT = 
ySol = iztrans(pZT,z,k) % se aplica la transformada inversa, la variable z pasa a k
ySol = 
ySol = simplify(ySol) % de ser posible se simplifica la expresión resultante, esta es la solución de la ecuación en diferencias
ySol = 
Gráfica de la solución analítica
kValues = 0:10; %instantes de muestreo
T=0.2;% intervalo de muestreo
ySolValues = subs(ySol,k,kValues); %evaluamos los valores de k en la solución
ySolValues = double(ySolValues); % pasamos los valores de simbólico a variable numérica
ySolValues = real(ySolValues); % capturamos la parte real de los valores
figure(1)
clf
stem(kValues*T,ySolValues) %gráfica de la función discretizada
title('Solution of equation')
xlabel('t')
ylabel('y(kT)')
grid on

Solución mediante la serie infinita de potencias

Para resolver la ecuación mediante la serie infinita de potencias debemos aplicar la transformada Z y llegar a la expresión racional . Para el caso del ejemplo que estamos resolviendo tenemos que . Podemos aplicar la división larga para obtener los primeros 11 valores de la función (equivalentes a los puntos de la gráfica anterior) mediante la función ldiv desarrollada por Tamer Melik. El primer argumento que recibe esta función son los coeficientes del numerados, luego los coeficientes del denominador y finalmente el número de valores de la función por obtener
yvals=ldiv([5,4,5],[5,4],11)
yvals = 1×11
1.0000 0 1.0000 -0.8000 0.6400 -0.5120 0.4096 -0.3277 0.2621 -0.2097 0.1678

Clases 7 y 8: Introducción a las funciones de transferencia

Definición de la función de transferencia

Cuando aplicamos la transformada de Laplace a un problema de valores iniciales no homogéneo
obtenemos una expresión de la forma
Adicionalmente, si la entrada fuera de la forma , llegaríamos a que
Donde es el polinomio característico de la ecuación diferencial homogénea en la variable s, es un polinomio generado a partir de las condiciones iniciales y es un polinomio generado a partir de las derivadas de la entrada (las condiciones iniciales de la entrada se asumen iguales a 0). De la última expresión que obtuvimos podemos observar que la solución de la ecuación diferencial no homogénea viene dada por una respuesta a las condiciones iniciales y una respuesta a la entrada . Es decir, la solución de la ecuación se compone de la suma de la solución a la ecuación asociada homogénea y una solución particular con condiciones iniciales nulas.
De la naturaleza de la transformada de Laplace se sigue que y que .
Para el caso discreto, de manera análoga, la solución mediante transformada Z de la ecuación en diferencias no homogénea con condiciones iniciales
nos lleva a la forma
Y de la naturaleza de la transformada Z se sigue que y que .
Se hace de particular relevancia, el caso en el que las condiciones iniciales son nulas puesto que la expresión para se torna de la forma
A la expresión se le conoce como función de transferencia del sistema de ecuaciones diferenciales / en diferencias. La función de transferencia es de particular relevancia para el área de sistemas dinámicos, donde podemos resignificar cada componente de la expresión anterior en términos de una planta/proceso/sistema (dado por ), un conjunto de señales de entrada que se aplican al sistema (), un conjunto de estados del sistema () y un conjunto de salidas del sistema (). El tema de los estados del sistema será abordado más adelante. Lo anterior, nos permite representar el anterior problema matemático de la forma
En este orden de ideas, la función de transferencia es una caraterización completa del sistema que nos permite determinar la salida del sistema para cualquier entrada que se le aplique. En términos matemáticos, esto puede verse al notar que . Es decir que si conocemos la función de transferencia del sistema y la entrada que se le aplica, podemos calcular la salida mediante transformada inversa o convolución.
En general, el concepto de función de transferencia también es aplicable a cualquier sistema en estado estacionario (sin variaciones) aunque las condiciones de dicho estado no sean nulas. Para ello, basta trabajar con las variables incrementales . Es decir, obtenemos el cambio que experimenta la salida del sistema al modificar la entrada en estado estacionario y luego adicionamos el valor de la salida en estado estacionario para obtener el comportamiento del sistema, tal y como se ejemplifica en el siguiente diagrama:

Conceptos claves

  1. Polos: Son los ceros del polinomio característico , por tanto, son los valores de para los cuales la función de transferencia tiende a infinito. Los polos determinan la estabilidad del sistema. Para el caso continuo, si todos los polos se encuentran en el semiplano izquierdo del plano s se sigue que el sistema será estable. Para el caso discreto, un sistema será estable si todos los polos se encuentran dentro del círculo unitario del plano z.
  2. Ceros finitos: Son los ceros del polinimio de la entrada , por tanto, son los valores de para los cuales la función de transferencia tiende a cero. Los ceros nos permiten identificar entradas que son anuladas por el sistema e intuir ciertos comportamientos del sistema, tal y como se verá más adelante.
  3. Orden: Es el grado de , dicho grado nos habla del orden de la ecuación diferencial asociada al sistema.
  4. Orden relativo: . El orden relativo nos habla sobre el retardo del sistema para el caso discreto. Este tema se abordará con más profundidad en la siguiente clase.
  5. Polos dominantes: Polos estables cercanos al eje imaginario (caso continuo). Polos estables cercanos al círculo unitario (caso continuo). Debido a dicha cercanía, estos polos son los que más determinan el comportamiento del sistema.
  6. Polos insignificantes: Polos alejados del polo más dominante. Se calcula que un polo ajelado más de 5 veces del polo dominante podría ser insignificante. Un polo alejado más de 10 veces del polo dominante se juzga insignificante. Los polos insignificantes aportan muy poco a la dinámica del sistema.
  7. Ecuación característica: .

Clases 9 y 10: Funciones de transferencia

Retardos y orden relativo

Las transformadas para una señal continua con retardo son de la forma
Por tanto, un problema como el siguiente
conduce a la solución
mientras que un problema como
conduce a la solución
Por tanto, la presencia de exponenciales, en el caso continuo, o la presencia de polos en , en el caso discreto, nos alertan sobre la presencia de retardos en el sistema. Para el caso continuo los retardos son bastante intuitivos mientras que el caso discreto requiere un poco más de atención. En primer lugar, ya no podremos interpretar inmediatamente que el denominador de la función de transferencia corresponde con el polinomio característico de la ecuación en diferencias asociada al sistema. Para determinar el polinomio característico tendremos que factorizar, de ser posible, el denominador, hasta que obtengamos una expresión de la forma
Ecuación en diferencias a partir de la función de transferencia
Partamos de un sistema con la función de transferencia , dado que cuando trabajamos con funciones de transferencia asumimos que las condiciones iniciales son iguales a 0, se sigue que la ecuación en diferencias debe ser de la forma . Ahora, hagamos y solucionemos este sistema empleando MatLab
clear;clc
G = tf([1 0],[1 0 0.25],1)
G = z ---------- z^2 + 0.25 Sample time: 1 seconds Discrete-time transfer function.
figure(1)
clf
step(G)
xlim([0 10])
ylim([0 1.5])
Sin embargo, debemos notar que MatLab toma un valor distinto de 0 para . Para entender porqué sucede esto, podemos modificar la función de transferencia para obtener la ecuación en diferencias hacia atrás
En este caso, podemos resolver , que coincide con la gráfica entregada por MatLab.
El orden relativo revisitado
Del ejemplo anterior podemos empezar a intuir el efecto del orden relativo, para el caso discreto, en el comportamiento del sistema. De hecho, el orden relativo de la función de transferencia es equivalente al retardo que presenta el sistema. Esto puede comprobarse cuando determinamos la ecuación en diferencias hacia atrás. Por ejemplo:
Notemos que el orden relativo de esta ecuación en diferencias es 0 y que , que coincide con la solución de MatLab
G = tf([1 0 0],[1 0 0.25],1)
G = z^2 ---------- z^2 + 0.25 Sample time: 1 seconds Discrete-time transfer function.
figure(2)
clf
step(G)
xlim([0 10])
ylim([0 1.5])

Sistemas de fase no mínima

El efecto de los ceros de la función de transferencia

Diremos que un sistema es de fase no mínima cuando tiene al menos un cero finito por fuera de la región de estabilidad en el plano S o Z, según sea el caso. Los sistemas de fase no mínima presentan comportamientos poco intuitivos debido al efecto de los ceros. En la sección anterior vimos cómo podíamos deducir la ecuación diferencial/diferencias a partir de la función de transferencia. Podemos notar que el componente determina la forma de la entrada para la ecuación diferencial/diferencias, es decir, el numerador de la función de transferencia afecta a la entrada del sistema. Pensemos, por ejemplo, en una función de transferencia dada por , la presencia de s en el numerador sugiere que la entrada que se aplique al sistema será derivada, ¿qué sucedería entonces si a este sistema la aplicamos una entrada constante ? Veámoslo mediante simulación
clear;clc;close all
G=zpk(0,[-1 -2 -3],1)
G = s ----------------- (s+1) (s+2) (s+3) Continuous-time zero/pole/gain model.
figure(1)
clf
step(G)
Vemos entonces que nuestro sistema tiene una reacción inicial al aplicar la entrada pero luego se estabiliza en 0, a pesar de que la entrada no es nula. ¿Qué ha pasado entonces? ¡La entrada ha sido anulada por el cero de la función de transferencia! Para analizar este proceso matemáticamente, primero veamos la forma de la ecuación diferencial que correspondería con este problema
tf(G)
ans = s ---------------------- s^3 + 6 s^2 + 11 s + 6 Continuous-time transfer function.
. Ahora bien, considerando que la entrada del sistema es el escalón, tenemos que y por tanto la ecuación correspondiente es , cuya solución corresponde con el valor de estabilización de la simulación. Sin embargo, se aprecia una dinámica al principio de la simulación que no corresponde con la la ecuación diferencial que hemos hallado. Esto se debe a que . Visto desde la perspectiva de las funciones de transferencia, el problema que estamos resolviendo toma la forma , resta notar que para que ambos resultados sean coherentes. Ahora bien, podemos continuar con la solución para haciendo fracciones parciales y aplicando transformada inversa
[r,p,k] = residue(1,[1 6 11 6])
r = 3×1
0.5000 -1.0000 0.5000
p = 3×1
-3.0000 -2.0000 -1.0000
k = []
Este resultado debe coincidir con la respuesta al impulso para la función de transferencia . Veamos
G2=tf(1,[1 6 11 6])
G2 = 1 ---------------------- s^3 + 6 s^2 + 11 s + 6 Continuous-time transfer function.
figure(2)
clf
subplot(1,2,1)
impulse(G2)
legend('G_2')
subplot(1,2,2)
step(G)
legend('G')
En conclusión, el cero de la función de transferencia actúa anulando la entrada que se le aplica al sistema pero dicha anulación no ocurre inmediatamente dado que en el instante inicial aparece un impulso. Lo particularmente interesante de este resultado es que los sistemas de fase no mínima pueden anular entradas que divergen asintóticamente, haciendo que el sistema siga convergiendo. En este orden de ideas, un sistema con los siguientes ceros anulará entradas de la forma .

Comportamiento de un sistema de fase no mínima

Además de las propiedades mencionadas anteriormente, los sistemas de fase no mínima tienen comportamientos interesantes. Si bien el nombre de fase no mínima corresponde a dinámicas que se observan en el análisis en frecuencia de estos modelos, una forma rápida de entender dicho nombre es pensando que su comportamiento es más rico en términos de la fase, es decir, del ángulo, es decir, más ricos en términos de dinámica. Algo que se observa con frecuencia en estos modelos, es que al aplicarles una entrada reaccionan primero en el sentido opuesto al esperado, aunque luego se comporten de forma esperada. Tal comportamiento hace que sea especialmente difícil diseñar un control para estos sistemas. Algunos ejemplos cotidianos de sistemas de fase no mínima son:
Veamos la respuesta a una entrada de tipo escalón para algunos sistemas de fase no mínima. Primero, analicemos un modelo con . En este caso, vemos que aunque estamos aplicando una entrada positiva, el sistema responde disminuyendo su valor en estado de equilibrio (esto se debe al efecto del numerador de la función de transferencia , que puede entenderse como derivar la entrada y luego restarle 3 veces el valor de la entrada. Notemos que el sistema muestra una ganancia al principio de la simulación.
clear;clc
G=zpk([3],[-1 -1.5 -2],1);
tf(G)
ans = s - 3 ------------------------- s^3 + 4.5 s^2 + 6.5 s + 3 Continuous-time transfer function.
figure(3)
clf
step(G)
Aumentemos ahora el número de ceros del sistema
G2=zpk([3 2],[-1 -1.5 -2],1);
tf(G2)
ans = s^2 - 5 s + 6 ------------------------- s^3 + 4.5 s^2 + 6.5 s + 3 Continuous-time transfer function.
figure(4)
clf
step(G2)
En este caso, el sistema presenta una oscilación al inicio de la dinámica. Podría pensarse que hay presencia de valores propios imaginarios pero no es el caso. El sistema se estabiliza ahora en un valor mayor al que tenía en estado estacionario (puede hacerse el mismo análisis del caso anterior para comprender por qué). Notamos entonces que la presencia de ceros de fase no mínima enriquece la dinámica del sistema. Podemos comprobar tal afirmación simulando un sistema con el mismo polinomio característico pero de fase mínima
G3=zpk([-3 -4],[-1 -1.5 -2],1);
tf(G3)
ans = s^2 + 7 s + 12 ------------------------- s^3 + 4.5 s^2 + 6.5 s + 3 Continuous-time transfer function.
figure(5)
clf
step(G3)

Reducción de orden

En general, es posible aproximar un modelo de orden superior con uno de menor orden que conserve sus principales características. Existen varios métodos formales de reducción de orden pero en este curso nos concentraremos en el método de cancelación de polos y ceros de fase mínima y el método de eliminación de polos insignificantes. Estos métodos explotan propiedades matemáticas intuitivas de la estructura de la función de transferencia para eliminar polos y luego compensar dicha eliminación con una ganancia , de forma que se conserve el valor de estabilización del sistema ante una entrada escalón. Para esta última parte, basta recordar que

Cancelación de polos y ceros de fase mínima

Pensemos que tenemos una función de transferencia con un polo y un cero cercanos
Para notar cuándo podemos decir que debemos ver que
.
De este resultado podemos concluir que el único caso en el que es posible cancelar el polo y el cero ocurre cuando el polo se encuentra en la región de estabilidad, es decir, . Por otro lado, la similitud entre y dependerá de la cercanía a 0 de la diferencia entre el polo y el cero y de velocidad con la que el término desaparece . El caso discreto es análogo al continuo, se concluye que el polo debe estar dentro de la región de estabilidad. Ahora bien, cancelar un polo y un cero cuando el polo es estable y el cero se encuentra fuera de la región de estabilidad, causaría que un sistema de fase no mínima se convirtiera en uno de fase mínima, lo cual no tiene mucho sentido.

Eliminación de polos insignificantes

De las secciones anteriores, sabemos que un polo es insignificante cuando es estable y está alejado de 5-10 veces o más del polo dominante (en norma). Pensemos entonces en una función de transferencia estable con un polo insignificante de la forma
, siendo la métrica inducida por la norma euclídea. Sabemos que es posible plantear una descomposición en fracciones parciales de la forma
. Ahora bien, la eliminación se basa en que la contribución de es despreciable con relación al resto de términos que aparecen en la descomposición de fracciones parciales. También es posible proponer que , siempre y cuando la fracción siga siendo propia.

Ejemplos de reducción de orden

Ejemplo 1
Supongamos que tenemos la siguiente función de transferencia
clear;clc
syms s
G=zpk([-3 -4],[-1 -10 -2],1)
G = (s+3) (s+4) ------------------ (s+1) (s+2) (s+10) Continuous-time zero/pole/gain model.
figure(1)
clf
step(G)
En primer lugar, podemos plantear que , donde , es decir que . Pero también podríamos emplear una construcción más elaborada calculando las fracciones parciales de
f=(s+3)*(s+4)/(s^3+13*s^2+32*s+20)
f = 
partfrac(f)
ans = 
En este punto, despreciamos la fracción paracial que queremos remover y reconstruimos la expresión con orden reducido
expre=simplify(partfrac(f)-7/(12*(s+10)))
expre = 
[Eqnn,Eqnd] = numden(expre)
Eqnn = 
Eqnd = 
Con esta información podemos plantear que , donde , es decir que
G2=tf(6/13*[1 13/5],[1 3 2]);
G1=zpk([-3 -4],[-1 -2],1/10);
figure(1)
clf
step(G)
hold on
step(G1)
step(G2)
legend({'G','G_1','G_2'})
Ejemplo 2
Supongamos que tenemos la siguiente función de transferencia
clear;clc
G=zpk([0.8 -0.4],[0.81 0.1-0.3i 0.1+0.3i],1,1)
G = (z-0.8) (z+0.4) --------------------------- (z-0.81) (z^2 - 0.2z + 0.1) Sample time: 1 seconds Discrete-time zero/pole/gain model.
figure(1)
clf
step(G)
Cancelando el cero y el polo de fase mínima cercanos tenemos que . Ahora, para calcular k hacemos
syms z k
f(z)=(z-0.8)*(z+0.4)/((z-0.81)*(z^2-0.2*z+0.1));
f1(z)=(z+0.4)/(z^2-0.2*z+0.1)*k;
solve(f(1)==f1(1),k)
ans = 
G1=tf(20/19*[1 0.4],[1 -0.2 0.1],1)
G1 = 1.053 z + 0.4211 ----------------- z^2 - 0.2 z + 0.1 Sample time: 1 seconds Discrete-time transfer function.
figure(2)
clf
step(G)
hold on
step(G1)
legend({'G','G_1'})

Caso especial de reducción de orden

Los métodos de reducción de orden que hemos visto no son aplicables cuando los ceros de la función de transferencia reducida actúan como anuladores de una entrada tipo escalón dado que no es posible identificar la ganancia. Lo máximo a lo que llegaremos en estos casos es proponer .

Función de ponderación (Secuencia de ponderación)

Pensemos que tenemos una función de transferencia estable y le aplicamos un impulso como entrada, . Por el teorema de convolución, la respuesta del sistema estará dada por , siendo . Ahora, dado que el sistema es estable, se sigue que y por ende, podemos aproximar esta secuencia haciendo . Llamaremos a como la secuencia de ponderación y a como la secuencia truncada de ponderación. Debe notarse que es una caracterización completa del sistema, es decir, si conocemos la secuencia de ponderación del sistema entonces podemos calcular la respuesta del sistema a cualquier entrada haciendo convolución. De forma tal, la secuencia truncada de ponderación es una aproximación al comportamiento del sistema. Es posible realizar este mismo análisis para sistemas de tiempo contínuo con la diferencia de que no obtendremos secuencias de ponderación sino funciones de ponderación. Las secuencias de ponderación son más versátiles puesto que es posible obtener los n términos de desde mediante división larga. En cualquier caso, tanto la secuencia como la función de ponderación, corresponden con la respuesta del sistema a una entrada de impulso, por tanto, se les conoce como respuesta al impulso.

Aproximación de la respuesta temporal mediante la secuencia de ponderación

Supongamos que tenemos la siguiente función de transferencia
clear;clc
G=zpk([1.01 -0.4],[0.81 0.1-0.3i 0.1+0.3i],1,1)
G = (z-1.01) (z+0.4) --------------------------- (z-0.81) (z^2 - 0.2z + 0.1) Sample time: 1 seconds Discrete-time zero/pole/gain model.
figure(1)
clf
step(G)
Para calcular los primeros 15 términos de la secuencia de ponderación truncada empleamos la función ldiv
tf(G)
ans = z^2 - 0.61 z - 0.404 -------------------------------- z^3 - 1.01 z^2 + 0.262 z - 0.081 Sample time: 1 seconds Discrete-time transfer function.
[h_n]=ldiv([0 1 -0.61 -0.404],[1 -1.01 0.262 -0.081],15)
h_n = 1×15
0 1.0000 0.4000 -0.2620 -0.2884 -0.1903 -0.1378 -0.1127 -0.0931 -0.0757 -0.0612 -0.0495 -0.0401 -0.0325 -0.0263
Con esta secuencia de ponderación podemos estimar fácilmente la respuesta del sistema al escalón unitario. Esto es particularmente simple dado que, por la estructura de la convolución, la respuesta al escalón unitario es la suma acumulada de la secuencia de ponderación
y_step=cumsum(h_n);
Ahora bien, dado que la secuencia de ponderación es la respuesta al impulso del sistema. Podemos comparar la secuencia de ponderación truncada con la secuencia de ponderación haciendo
figure(2)
clf
impulse(G)
hold on
stem(0:14,h_n)
legend({'G','h_n'})
Finalmente, comparamos la respuesta estimada al escalón unitario con la respuesta real del sistema haciendo
figure(3)
clf
step(G)
hold on
stem(0:14,y_step)
legend({'G','G estimada'})

Discretización de la función de transferencia

Transformada z modificada

La transformada z modificada (o avanzada) es una extensión de la transformada z que nos permite incorporar retardos que no son múltiplos del periodo de muestreo. La transformada z modificada tiene la forma
siendo
Si consideramos que el retardo es fijo, entonces todas las propiedades de la transformada z aplican para la transformada z modificada.

Ejemplos de discretización

Supongamos que tenemos un modelo dado por la siguiente función de transferencia
clear;clc
G=tf([1],[6 1],'InputDelay',0.45)
G = 1 exp(-0.45*s) * ------- 6 s + 1 Continuous-time transfer function.
Ahora, si queremos discretizar este modelo con un un petriodo de muestreo tal que el retardo no sea múltiplo, por ejemplo , entonces tendremos que emplear la transformada z modificada. Por suerte, MatLab puede realizar fácilmente el proceso de discretización mediante la función c2d.
Gd=c2d(G,0.2)
Gd = 0.02469 z + 0.008094 z^(-3) * -------------------- z - 0.9672 Sample time: 0.2 seconds Discrete-time transfer function.

Clase 11: Introducción al espacio de estado

Definiciones

  1. Estado: Es una caracterización completa del sistema. Consta de las ecuaciones de estado y las variables de estado. A veces, utilizaremos la palabra estado para referirnos al valor de las variables de estado en determinado instante de tiempo (también podemos llamar a esto el valor de estado).
  2. Variables de estado: Conjunto mínimo de variables necesarias para describir el estado.
  3. Espacio de estado: Espacio dimensional que contiene todos los posibles valores de estado del sistema. Una representación natural para el espacio de estado se obtiene al escoger como ejes de coordenadas a las variables de estado.
  4. Ecuación de estado: Sistema de n ecuaciones diferenciales de primer orden que involucran a las n variables de estado. La ecuación de estado determina un campo vectorial al interior del espacio de estado, de forma tal, conocer el valor de estado para determinado momento junto con las ecuaciones de estado permite conocer el comportamiento pasado y futuro del sistema.
  5. Vector de estado: Vector formado por las variables de estado. .
  6. Trayectoria de estado: Colección de todos los valores asumidos por el vector de estado para un intervalo de tiempo. Puede entenderse como una curva paramétrica (trayectoria) en el espacio de estado.

Ecuaciones de estado

Ecuación en el espacio de estado continuo

Una ecuación de espacio de estado continuo tiene la forma
con , además, definiendo los vectores
Podemos reescribir la ecuación de estado de la forma

Ecuación en el espacio de estado discreto

De manera similar al espacio de estado continuo, las ecuaciones en el espacio de estado discreto se pueden escribir de la forma
siendo
Cuando las ecuaciones de estado son lineales las funciones f y g se pueden desagregar en funciones lineales para escribirlas de la forma
o , aunque normalmente, (el sistema es propio).

Obtención de la ecuación de estado a partir de la ecuación diferencial

Es posible llevar una ecuación de orden n a una ecuación de estado mediante un cambio de variables.

Para sistemas continuos

Partiendo de la ecuación diferencial
Definimos n variables de estado de la forma
Por tanto, si la ecuación original es lineal puede reescribirse como un sistema n ecuaciones lineales de primer orden de la forma
Resta notar que normalmente no estamos interesados en conocer el comportamiento de todas las derivadas en una ecuación diferencial. En general, encontraríamos soluciones de la forma , por tanto, la salida de nuestro sistema correspondería con , es decir que
Caso especial: Presencia de derivadas en la entrada
Si partimos de una ecuación de la forma
Entonces tendríamos problemas para representar el sistema en la forma . Por tal razón, para obtener la ecuación de estado en este caso se añaden derivadas de la entrada en las ecuaciones de estado, empezando por la última ecuación de estado, de forma que se anulen las derivadas indeseadas. Los coeficientes de los términos que se añaden deben ser determinados cuando se calculan las derivadas de las variables de estado. En general, tendremos que añadir términos para eliminar las derivadas de u. Veamos los siguientes ejemplos.
1)
Para esta ecuación diferencial, la ecuación de estado resultante toma la forma
2)
Para esta ecuación diferencial, la ecuación de estado resultante toma la forma

Para sistemas discretos

Se procede de manera análoga al caso de los sistemas continuos. Partiendo de la ecuación en diferencias
se definen n variables de estado de la forma
Por tanto, si la ecuación original es lineal puede reescribirse como un sistema n ecuaciones lineales de primer orden de la forma
.
Caso especial: Presencia de diferencias en la entrada
Si partimos de una ecuación de la forma
Entonces añadiremos términos, de forma apropiada, para eliminar las diferencias de u. Para ver esto a través de un ejemplo, supongamos que queremos llevar a la representación en espacio de estado a la siguiente ecuación en diferencias
Primero, proponemos las variables de estado
De la ecuación para resulta inmediato que , luego, de la ecuación para concluimos primero que y con esto concluimos que . De forma tal
y así
Nota: En caso de que la ecuación en diferencias esté dada hacia atrás, entonces se debe hacer el cambio de variables para que quede en diferencias hacia adelante

Solución de la ecuación lineal de estado continua

En esta sección veremos métodos de solución para el problema de la forma

Solución mediante valores y vectores propios

Nos proponemos buscar una solución para el problema homogéneo de la forma , de existir tal solución, tendremos entonces que
, es decir que K es un autovector de A asociado al valor propio λ. Ahora bien, cuando solucionamos este problema de autovalores, podemos encontrarnos con 3 casos.
  1. Autovalores con multiplicidad algebráica igual a su multiplicidad geométrica. En este caso, es posible encontrar tantos autovectores linealmente independientes como dimensión tenga A. Luego, basta recordar que una matriz cuyas columnas están formadas por las soluciones asociadas a estos autovectores es una matriz fundamental para el sistema de ecuaciones diferenciales.
  2. Autovalores complejos. Los autovalores complejos siempre vienen en pares conjugados. Se tiene, además, que si es autovector asociado al autovalor complejo , entonces es autovector asociado con el autovalor complejo . Luego, recordando que la combinación lineal de soluciones de la ecuación es también solución de la ecuación, planteamos dos soluciones linealmente independientes de la forma .
  3. Autovalores con multiplicidad geométrica menor a la multiplicidad algebráica. En este caso, debemos resolver un problema de vectores propios generalizados. Cuando tenemos un autovector que corresponde con cierto autovalor para el cual queremos obtener m soluciones linealmente independientes es posible demostrar que dichas soluciones tienen la forma , donde son una cadena de vectores propios generalizados y son las m soluciones linealmente independientes. Ahora bien, la cadena de vectores propios generalizados se obtiene resolviendo .
Una vez se obtiene el conjunto con las m soluciones linealmente independientes se plantea la matriz fundamental y se resuelve el problema no homogéneo por el método de variación de parámetros:

Solución mediante serie de potencias (matriz exponencial o de transición del estado)

Definimos la matriz exponencial asociada a A como . En clase se demostró que es matriz fundamental principal para el sistema de ecuaciones diferenciales. De forma tal, una vez se le conoce es posible utilizar la versión simplificada del método de variación de parámetros:
A la matriz exponencial también se le conoce también como matriz de transición del estado y cuenta con las siguientes propiedades (Que se demostraron en clase). Sea
Obtener la matriz exponencial no suele ser sencillo salvo algunos casos. Por ejemplo, si la matriz A es diagonal, se sigue que
.
Otro caso relativamente sencillo ocurre cuando la matriz A es nilpotente, es decir que . En este caso, puede obtenerse como una sumatoria de finitos términos (que aún pueden ser muchos). Es sabido que todos los autovalores de una matriz nilpotente son 0, es decir que la ecuación característica de una matriz nilpotente es de la forma , por tanto, a lo sumo tendremos que calcular términos de la sumatoria. Por ejemplo, todos los autovalores de la matriz
clear;clc
A=sym([0 2 1;0 0 4; 0 0 0])
A = 
son 0
eig(A)
ans = 
Por tanto la matriz es nilpotente:
A^2
ans = 
Y así,
Un caso que es similar al de la matriz nilpotente es cuando el polinomio característico de la matriz consta de un único autovalor . En este caso podemos explotar algunas propiedades de la matriz de transición del estado:
Esta última expresión es de utilidad puesto que, por el teorema de Cayley-Hamilton, la matriz A satisface el polinomio característico, es decir . Por tanto, la matriz es nilpotente. Veamos el siguiente ejemplo
A=sym([2 1 1; 1 2 1; -2 -2 -1])
A = 
eig(A)
ans = 
Como la matriz tiene un único autovalor, proponemos
y comprobamos el comportamiento de la matriz
M=A-eye(3)
M = 
M^2
ans = 
De forma tal
De una manera más general, para obtener la matriz de transición del estado tendremos que diagonalizar la matriz o llevarla a la forma canónica de Jordan. Para ello, debemos recordar que una matriz es diagonalizable siempre que todos sus autovalores tengan multiplicidad algebráica igual a su multiplicidad geométrica. En caso tal, la matriz formada por los autovectores de la matriz que queremos diagonalizar hará el proceso de diagonalización. Veamos el siguiente ejemplo:
Dada la siguiente matriz
A=sym([2 1 1; 1 2 1; 2 -2 -1])
A = 
Calculamos los valores y vectores propios haciendo
[vecs,vals]=eig(A)
vecs = 
vals = 
Y diagonalizamos la matriz haciendo
inv(vecs)*A*vecs
ans = 
La diagonalización de la matriz A es particularmente útil puesto que, como se puede demostrar fácilmente, y por tal razón . Así, en el ejemplo que estamos realizando, obtenemos que
Podemos comprobar este resultado en MatLab haciendo
D=sym(diag([exp(3*t) exp(t) exp(-t)]))
D = 
vecs*D*inv(vecs)
ans = 
O directamente calculando la matriz exponencial de A con el comando expm.
expm(A*t)
ans = 
En el caso más general, resulta que cualquier matriz puede llevarse a la forma canónica de Jordan. Una matriz en forma de Jordan tiene sus valores propios en la diagonal y valores de 1 sobre la diagonal para aquellos valores propios cuya multiplicidad algebráica sea menor a la geométrica. La matriz que nos permite realizar esta transformación está formada por los vectores propios generalizados, los cuales, a su vez, se obtienen al resolver una cadena de la forma
Tomemos un ejemplo que ya resolvimos anteriormente
A=sym([2 1 1; 1 2 1; -2 -2 -1])
A = 
y veamos los valores y vectores propios
[vecs,vals]=eig(A)
vecs = 
vals = 
De aquí podemos notar que la multiplicidad algebráica es mayor a la geométrica. Exploremos ahora el comportamiento de
M=A-eye(3)
M = 
rank(M)
ans = 1
M^2
ans = 
Con lo anterior notamos que el comportamiento de M es problemático. Debemos encontrar 1 vector linealmente independiente adicional pero el rango y el índice de nulidad nos traerán problemas. Si quisiéramos resolver el problema mediante la cadena de vectores generalizados, podríamos proponer que
, pero este problema no tiene solución ya que el rango de M es el espacio generado por el vector . Ahora bien, dado que , podemos iniciar la cadena de vectores propios generalizados partiendo desde cualquier vector. Por ejemplo
, con lo cual y son vectores propios generalizados. Para encontrar el tercer vector propio, nos dirigimos a los resultados del problema de autovalores original y elegimos o , de forma tal, concluímos que las matrices o nos permiten llevar a la matriz A a la forma de Jordan. Veamos
Q1=sym([1 1 -1;1 0 1;-2 0 0]);
Q2=sym([1 1 -1;1 0 0;-2 0 1]);
inv(Q1)*A*Q1
ans = 
inv(Q2)*A*Q2
ans = 
Con eso, hemos obtenido una matriz Q tal que . La utilidad de encontrar la forma de Jordan para una matriz radica en que J puede escribirse como la suma entre una matriz diagonal y una nilpotente , donde, adicionalmente, D y N conmutan. Así, explotando la propiedad (6) de la matriz de transición del estado, tenemos que
Y por tanto
Que coincide con la respuesta que se halló con anterioridad. Este resultado puede comprobarse con MatLab haciendo
Q1*[1 t 0;0 1 0;0 0 1]*inv(Q1)*exp(t)
ans = 

Solución mediante transformada de Laplace

Es posible aplicar la transformada de Laplace para obtener la solución general del problema rápidamente. Veamos
Ahora, decimos que para obtener
Pero entonces debe ser matriz fundamental principal (de transición del estado) y por propiedades (2) y (3) se llega a que:

Ejemplo

clear;clc
Supongamos que queremos solucionar una ecuación de estado de la forma
Primero calcularemos la matriz de transición del estado.
syms t s c_1 c_2 c_3
A=sym([3 0 1;1 2 -1;-1 0 3])
A = 
Para decidir qué método elegir, calculamos los valores y vectores propios
[vecs,vals]=eig(A)
vecs = 
vals = 
Notamos que tenemos valores propios imaginarios. Es posible calcular la matriz exponencial pero luego tendríamos que extraer la parte real, lo cual puede ser tedioso. Para evitar este proceso, decidimos seguir el procedimiento indicado para resolver la ecuación homogénea por el método de valores y vectores propios. Encontraremos entonces 3 soluciones linealmente independientes para el espacio de estado. El valor propio real nos arroja la solución
x_1=sym([0 1 0]'*exp(2*t))
x_1 = 
Para obtener las soluciones linealmente independientes asociadas a los valores propios complejos procedemos haciendo
K=sym([1i -1 1]')
K = 
B_1=real(K)
B_1 = 
B_2=imag(K)
B_2 = 
x_2=(B_1*cos(t)-B_2*sin(t))*exp(3*t)
x_2 = 
x_3=(B_2*cos(t)+B_1*sin(t))*exp(3*t)
x_3 = 
De tal forma, la matriz fundamental estará dada por
phi(t)=[x_1 x_2 x_3]
phi(t) = 
Y la matriz fundamental principal (de transición del estado) estará dada por
psi(t)=phi(t)*inv(phi(0))
psi(t) = 
Para solucionar el problema no homogéneo aplicamos el método de variación de parámetros
invpsi(t)=inv(psi(t))
invpsi(t) = 
x(t)=psi(t)*[1 -1 2]'+psi(t)*int(invpsi(t)*([1 0 1]'*sin(3*t)),[0 t])
x(t) = 
Resolviendo llegamos a que
x(t)=simplify(x(t))
x(t) = 
Ahora bien, también podríamos solucionar este problema mediante transformada de Laplace. Para ello hacemos
M=inv(s*eye(3)-A)
M = 
psi(t)=ilaplace(M,t)
psi(t) = 
Nótese que esta matriz de transición del estado es la misma que la obtenida por el método de valores y vectores propios. Aplicando el resto del método llegamos a que
invpsi(t)=inv(psi(t));
x(t)=simplify(psi(t)*[1 -1 2]'+psi(t)*int(invpsi(t)*([1 0 1]'*sin(3*t)),[0 t]))
x(t) = 
Para graficar es recomendable crear una función
sol=matlabFunction(x)
sol = function_handle with value:
@(t)[cos(t.*3.0).*(-3.3e+1./3.25e+2)-sin(t.*3.0).*(5.6e+1./3.25e+2)+exp(t.*3.0).*cos(t).*(3.58e+2./3.25e+2)+exp(t.*3.0).*sin(t).*(7.19e+2./3.25e+2);cos(t.*3.0).*(-6.0./3.25e+2)+exp(t.*2.0).*(1.6e+1./1.3e+1)+sin(t.*3.0).*(8.0./3.25e+2)-exp(t.*3.0).*cos(t).*(7.19e+2./3.25e+2)+exp(t.*3.0).*sin(t).*(3.58e+2./3.25e+2);cos(t.*3.0).*(-6.9e+1./3.25e+2)-sin(t.*3.0).*(5.8e+1./3.25e+2)+exp(t.*3.0).*cos(t).*(7.19e+2./3.25e+2)-exp(t.*3.0).*sin(t).*(3.58e+2./3.25e+2)]
t=0:0.001:2;
Y=eye(3)*sol(t);%y(t)=Cx(t)
figure(1)
clf
subplot(3,1,1)
plot(t,Y(1,:))
title('y_1')
subplot(3,1,2)
plot(t,Y(2,:))
title('y_2')
subplot(3,1,3)
plot(t,Y(3,:))
title('y_3')
Comprobación de resultados
Para comprobar los resultados empleamos las funciones ss y lsim de MatLab de la siguiente forma
1) Definimos las matrices del sistema en espacio de estado
A=[3 0 1;1 2 -1;-1 0 3];
B=[1 0 1]';
C=eye(3);
D=[0 0 0]';
x0=[1 -1 2]';
2) Creamos el modelo de espacio de estado
S=ss(A,B,C,D)
S = A = x1 x2 x3 x1 3 0 1 x2 1 2 -1 x3 -1 0 3 B = u1 x1 1 x2 0 x3 1 C = x1 x2 x3 y1 1 0 0 y2 0 1 0 y3 0 0 1 D = u1 y1 0 y2 0 y3 0 Continuous-time state-space model.
3) Creamos la entrada que aplicaremos al sistema
t=0:0.001:2;
u=sin(3*t);
4) Simulamos la respuesta del sistema a la respectiva entrada
figure(2)
clf
lsim(S,u,t,x0)
out=lsim(S,u,t,x0);
Gráfica para la primera variable de estado del modelo. Puede apreciarse que son idénticas.
figure(3)
clf
plot(t,out(:,1))
hold on
plot(t,resp(1,:))

Solución de la ecuación lineal de estado discreta

La ecuación de estado discreta se comporta de manera considerablemente distinta a la ecuación de estado continua. En general, excepto que la matriz A sea diagonalizable, será bastante complicado encontrar una solución mediante valores y vectores propios. Por otro lado, la naturaleza de la ecuación de estado discreta sugiere en sí misma un mecanismo iterativo para resolverse considerando que A es una matriz de transición del estado. Tomemos la forma general del espacio de estado discreto
y veamos cómo podemos resolverla

Solución iterativa

Podemos aprovecharnos de la estructura del espacio de estado discreto para notar que
Ahora, haciendo y sustituyendo la variable n resultante por la variable k, nuevamente, obtenemos la fórmula
Esta expresión, sin embargo, sigue siendo difícil de resolver. Deberíamos enfocar nuestros esfuerzos entonces en obtener una expresión analítica que nos de cuenta del comportamiento de . Podemos calcular dicha expresión aplicando la transformada Z. Veamos.

Solución mediante transformada Z

Aplicamos la transformada Z sobre para obtener
Ahora, de forma similar a lo que hicimos para la transformada de Laplace, definimos y así, tras aplicar la transformada Z inversa sobre la ecuación, llegamos a que
Claramente, aquí hemos obtenido una matriz que tiene ventajas con respecto al cálculo de la matriz , sin embargo, la convolución puede seguir siendo problemática. Ventajosamente, no es necesario realizar la convolución cuando se pueda calcular directamente a ya que la ecuación tendría la forma

Ejemplo

Resolvamos la siguiente ecuación de estado discreta
clear;clc
syms z k
Primero definimos las matrices involucradas, las condiciones iniciales y la entrada
A=[-0.1 1;0 0.2];
B=[1 1]';
C=eye(2);
D=[0 0]';
u=1;
x0=[1 -1]';
Luego definiremos las matrices a las que debemos calcular la transformada Z inversa
M=inv(z*eye(2)-A)*z
M = 
N=M/z*B*ztrans(u,z)
N = 
Y procederemos a calcular la transformada Z inversa para obtener y
phi(k)=iztrans(M,k)*x0
phi(k) = 
psi(k)=iztrans(N,k)
psi(k) = 
Finalmente, obtenemos
x(k)=simplify(phi(k)+psi(k))
x(k) = 
Para graficar es recomendable crear una función
sol=matlabFunction(x)
sol = function_handle with value:
@(k)[(1.0./5.0).^k.*(-1.5e+1./2.0)+(-1.0./1.0e+1).^k.*(7.1e+1./1.1e+1)+4.5e+1./2.2e+1;5.0.^(-k).*(-9.0./4.0)+5.0./4.0]
k=0:6;
Y=C*sol(k);%y(k)=Cx(k)
figure(1)
clf
subplot(2,1,1)
stairs(k,Y(1,:))
title('y_1')
subplot(2,1,2)
stairs(k,Y(2,:))
title('y_2')
Comprobación de los resultados
Para comprobar los resultados empleamos las funciones ss y lsim de MatLab de la siguiente forma:
Sd=ss(A,B,C,D,1)%el último parámetro es el tiempo de muestreo
Sd = A = x1 x2 x1 -0.1 1 x2 0 0.2 B = u1 x1 1 x2 1 C = x1 x2 y1 1 0 y2 0 1 D = u1 y1 0 y2 0 Sample time: 1 seconds Discrete-time state-space model.
figure(2)
clf
lsim(Sd,ones(7,1),k,x0)

Transformaciones lineales y formas canónicas

De las secciones anteriores ya hemos visto la utilidad de las transformaciones lineales para encontrar soluciones a las ecuaciones de estado. Dichas transformaciones pueden emplearse para realizar un cambio de base para el espacio de estado. Recordemos que la base del espacio de estado está dada por las variables de estado, de forma tal, partiendo de podemos plantear una transformación conveniente sobre A de forma tal que la relación entrada-salida de la ecuación de estado se conserve. De manera general, pensemos que que , donde J es la forma de Jordan de la matriz A. Ya hemos visto que la forma canónica de Jordan es particularmente útil a la hora de aplicar los métodos de solución que hemos visto. Ahora solo nos resta reescribir el espacio de estado para que se conserve la relación entrada-salida
Este mismo esquema puede seguirse con cualquier transformación lineal biyectiva (que tenga inversa). Las transformaciones lineales para el espacio de estado continuo son exactamente iguales a las del caso discreto.

Ejemplo

clear;clc
A=[2 1 1; 1 2 1; 2 -2 -1]
A = 3×3
2 1 1 1 2 1 2 -2 -1
B=[1 0 1]';
C=eye(3);
D=[0 0 0]';
x0=[1 -1 2]';
Definimos el sistema en espacio de estado
S=ss(A,B,C,D);
Ahora, lo transformaremos a la forma de Jordan, para ello usamos la función ss2ss una vez identifiquemos a la matriz Q. El argumento que debemos pasar a esta función para indicar la transformación a realizar es .
[vecs,vals]=jordan(sym(A))
vecs = 
vals = 
ivecs=inv(vecs)
ivecs = 
Finalmente, transformamos el modelo haciendo
S2=ss2ss(S,double(ivecs))
S2 = A = x1 x2 x3 x1 3 0 0 x2 0 1 0 x3 0 0 -1 B = u1 x1 1 x2 1 x3 0 C = x1 x2 x3 y1 1 0 -0.25 y2 1 -1 -0.25 y3 0 1 1 D = u1 y1 0 y2 0 y3 0 Continuous-time state-space model.

Función de transferencia MIMO

En los ejemplos que hemos visto hasta ahora de espacio de estado hemos lidiado con matrices B sencillas que corresponden con una única entrada al modelo, no obstante, podemos definir B como una matriz de siendo n el número de estados del sistema y un número de entradas. En este caso, dejaría de ser un escalar y pasaría a ser un vector de entradas arbitrarias. Para el caso de múltiples salidas ya hemos presenciado la forma que toma , siendo p el número de salidas del sistema. Cuando un sistema tiene una única entrada y una única salida, se le conoce como SISO (single input - single output) mientras que cuando consta de múltiples entradas y salidas se le conoce como MIMO (multiple inputs - multiple outputs). Para el caso de sistemas MIMO es igualmente posible plantear una función de transferencia, pero en este caso, adquiere forma matricial, tal y como veremos a continuación.
Tenemos que la función de transferencia se define como la siguiente relación , de forma tal aplicando la transformada de Laplace (Z) a la ecuación de estado con condiciones iniciales nulas vemos que
De forma tal, los polos del sistema estarán dados por los valores de para los cuales . Esto ocurre cuando la matriz no es invertible, es decir, su determinante es 0. Por tanto, los polos del sistema corresponden con los valores propios de la matriz A. Los ceros del sistema corresponden con quellos valores para los cuales , es decir, corresponden con el espacio nulo de G. Con lo que hemos visto hasta ahora, podemos conncluir que los conceptos de ecuación diferencial (en diferencias), ecuación de estado y función de transferencia están relacionados entre sí. De hecho, ya hemos visto cómo podemos obtener una de estas caracterizaciones del sistema a partir de la otra. El caso que no hemos visto es cómo podemos obtener la ecuación de estado a partir de la función de transferencia. Este último procedimiento se hace utilizando gráficos de flujo de señal y la regla de Mason, pero no exploraremos estas herramientas por el momento.
MatLab nos permite pasar de la función de transferencia a la ecuación de estado y viceversa con las funciones ss y tf, tal y como veremos a continuación

Ejemplos

Obtener la función de transferencia a partir de la ecuación de estado
clear;clc
A=[-5 -1;3 -1];
B=[2 5; 1 0]';
C=[1 2;0 1];
D=[0]';
S=ss(A,B,C,D)
S = A = x1 x2 x1 -5 -1 x2 3 -1 B = u1 u2 x1 2 1 x2 5 0 C = x1 x2 y1 1 2 y2 0 1 D = u1 u2 y1 0 0 y2 0 0 Continuous-time state-space model.
G=tf(S)
G = From input 1 to output... 12 s + 59 1: ------------- s^2 + 6 s + 8 5 s + 31 2: ------------- s^2 + 6 s + 8 From input 2 to output... s + 7 1: ------------- s^2 + 6 s + 8 3 2: ------------- s^2 + 6 s + 8 Continuous-time transfer function.
Para obtener los ceros y polos de la función de transferencia hacemos
pole(G)
ans = 4×1
-4 -2 -4 -2
tzero(G)
ans = 2×1
-4.0000 -2.0000
eig(A)
ans = 2×1
-4 -2
Obtener la ecuación de estado a partir de la función de transferencia
G=tf([12 59],[1 6 8])
G = 12 s + 59 ------------- s^2 + 6 s + 8 Continuous-time transfer function.
S=ss(G)
S = A = x1 x2 x1 -6 -2 x2 4 0 B = u1 x1 4 x2 0 C = x1 x2 y1 3 3.688 D = u1 y1 0 Continuous-time state-space model.

Discretización en el espacio de estado

Nota: Para discretizar con MatLab empleamos la función c2d.

Modelado de sistemas no lineales

Características de los sistemas no lineales

Principio de superposición

Los sistemas no lineales no cumplen con el principio de superposición. Esto nos brinda una herramienta para comprobar, en la práctica, si una planta se comporta de forma lineal o no: modificar la entrada de la planta y registrar los valores en los que la salida se estabiliza; si los valores de estabilización comparados con los valores de la entrada no forman una línea recta, se concluye que el sistema es no lineal.

Puntos de equlibrio

Para encontrar los puntos de equilibrio de un sistema lineal resolvemos una ecuación de la forma para el caso continuo, o de la forma para el caso discreto. En cualquiera de los dos casos, si la matriz A es invertible, el problema solo tendrá una solución. Por otro lado, si la matriz no es invertible, el problema tendrá un continuo de soluciones (considerando que siempre existe la solución ). Contrario a lo anterior, los sistemas no lineales pueden presentar múltiples puntos de equilibrio cuya estabilidad es de carácter local.

Bifurcación

Por bifurcaciones nos referiremos a cambios en la estructura topológica del espacio de estado ocasionados por variaciones en los valores de los parámetros. Las bifurcaciones pueden alterar el número de puntos de equilibrio y su estabilidad.

Caos

Los sistemas no lineales de orden mayor a 3 pueden presentar caos. No existe una definición precisa de caos pero se considera que éste comple con tres condiciones:
  1. Sensibilidad a las condiciones iniciales: ante variaciones arbitrariamente pequeñas en las condiciones iniciales de los parámetros el sistema reacciona de forma impredeciblemente diferente.
  2. Transitividad topológica: sea la función que nos da la trayectoria del sistema caótico partiendo de las condiciones iniciales dadas por y sea B el atractor del sistema. La propiedad de transitividad topológica indica que .
  3. No perioricidad: Un sistema caótico puede tener órbitas periódicas dentro de sí pero éstas órbitas no son atractores, de forma tal, toda solución que no parta de una órbita no será periódica.

Respuesta a entradas sinusoidales

Los sistemas lineales responden a las entradas sinusoidales con salidas sinusoidales donde se altera la magnitud y la fase. Un sistema no lineal puede responder a una entrada sinusoidal con cualquier estructura de salida.

Comportamiento asintótico

Los sistemas no lineales pueden tender a infinito en periodos de tiempo finitos mientras que los sistemas no lineales pueden hacerlo solo para periodos de tiempo infinitos.

Funciones adicionales para el curso

function funresult=ldiv(a,b,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function (ldiv) : calculate inverse Z-transform by long division
% Author : Tamer mohamed samy abdelazim Mellik
% Contact information :
%Department of Electrical & Computer Engineering,
%University of Calgary,
%2500 University Drive N.W. ,
%Calgary, AB T2N 1N4 ,
%Canada .
% email :abdelasi@enel.ucalgary.ca
% email : tabdelaz@ucalgary.ca
% Webpage : http://www.enel.ucalgary.ca/~abdelasi/
% Date : 2-5-2002
% Version : 1.0.0
%Example
% This function like deconv but it help if the numerator less or equal degree of denominator
% if you have this function (It must arranged in terms of minus power of Z):
% 1
% G(z)= -----------------
% -1 -2
% ( 5 - Z - 3 Z )
% and you want to calculate long division or inverse Z transform :
% The numerator is a=[1] and the denominator is b= [5 -1 -3 ]
% call the function ldiv(a,b) to get the funresult 20 items (default)
% another example :
% -2 -3
% ( 5 - 3 Z + 4 Z )
% G(z)-----------------
% -1 -2
% (5 - Z - 3 Z )
% a=[5 0 -3 4] , b= [5 -1 -3 ] and you want the funresult 100 terms !
% ldiv(a,b,100)
% Note : The author doesn't have any responsibility for any harm caused by the use of this file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%a numerator
%b denominator
%default order of the filter == 20
funresult=[];
if nargin < 3
if nargin > 1
N=20;
else
disp('Usage: M = ldiv(a,b,N)')
disp('a:numerator , b denominator and N is the order of the resultant filter')
return
end
end
if size(a) < 1
disp('Error: numerator must at least have one element not empty')
return
end
if size(b) < 1
disp('Error: denominator must at least have one element not empty')
return
end
if b(1)==0
disp('Error: The first element of denominator must have nonzero value')
return
end
if size(b) < 2
funresult=a./b;
for i =length(funresult)+1:N
funresult(i)=0;
end
return
end
for i = length(a)+1:N
a(i)=0;
end
for i = 1 : N
funresult(i)=a(1)/b(1);
if length(a)>1
for k= 2:length(b)
if k > length(a)
a(k)=0;
end
a(k)=a(k)-funresult(length(funresult))*b(k);
end
for i = 1:length(a)-1
a(i)=a(i+1);
end
a=a(1:length(a)-1);
end
end
end